# Project Info---------------
#  Project: Outgroup Avoidance
#  Purpose: Analyze First Survey
#  Outputs: Figure 7-8, S5-6 and Table S9


library("tidyverse")
library("estimatr")
library("stdidx")
library("ggthemes")

#Read data-----
stdy2_dta <- read.csv("study2.csv")

#Main Text----
##Figure 7-----
###Panel A------
stdy2_dta %>% 
  ggplot(., aes(x = choice_pal_fac))+
  geom_bar(stat="count",
                 fill = "dodgerblue4")+
  facet_wrap(~options_fac)+
  labs(x = "Select to Read about Palestinians in Choice Task",
       y = "# of Respondents")+
  theme_bw()+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 18),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())


###Panel B------
ideology <-
  lm_robust(choice_pal~ stdize(left_right_id),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Ideology")

religiosity <-
  lm_robust(choice_pal~ stdize(as.numeric(rlg)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Religiosity")


age_group <-
  lm_robust(choice_pal~ stdize(as.numeric(age)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Age Group")


therm <-
  lm_robust(choice_pal~ stdize(as.numeric(arab_therm)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Palastenian \nThermometer")


sm_usage <-
  lm_robust(choice_pal~ stdize(as.numeric(sm_usage)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Frequency on \nSocial Media")


sm_expose <-
  lm_robust(choice_pal~ stdize(as.numeric(sm_exposure)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Value \nDiversity \nOnline")


speak_arabic <-
  lm_robust(choice_pal~ stdize(as.numeric(speak_arabic)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Arabic \nSpeaking")



edu <-
  lm_robust(choice_pal~ stdize(as.numeric(education)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Education")

male <-
  lm_robust(choice_pal~ stdize(as.numeric(male)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Male")


covs <-
  rbind(ideology,religiosity,age_group,therm,sm_usage,sm_expose,speak_arabic,
        edu, male)


ggplot(covs, aes(x=estimate, y=reorder(term, estimate))) + 
  geom_point(color = "dodgerblue4")+
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high), 
                width=.1, color = "dodgerblue4") +
  geom_vline(xintercept = 0,
             color = "black",
             linetype = "dashed",
             size = 0.4) +
  xlim(-.2,.2)+
  labs(x = "Bivariate Correlation",
       y = "")+
  scale_color_manual(values = c("dodgerblue4"))+
  theme_bw()+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 16),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())


##Figure 8------
# reshape data
why_avoid <-
  as.data.frame(
    rbind(
      cbind(sum(grepl(1, stdy2_dta$post_task2)),"Feel \nUncomfortable"),
      cbind(sum(grepl(2, stdy2_dta$post_task2)),"Feel \nAnger"),
      cbind(sum(grepl(3, stdy2_dta$post_task2)),"Feel \nShame"),
      cbind(sum(grepl(4, stdy2_dta$post_task2)),"Feel \nDespair"),
      cbind(sum(grepl(5, stdy2_dta$post_task2)),"Reading \nwill \nnot \nChange \nAnything"),
      cbind(sum(grepl(6, stdy2_dta$post_task2)),"I've Seen \nEnough"),
      cbind(sum(grepl(7, stdy2_dta$post_task2)),"Other \nPriorities"),
      cbind(sum(grepl(8, stdy2_dta$post_task2)),"Other \nReason"),
      cbind(sum(grepl(9, stdy2_dta$post_task2)),"Feel \nFear"),
      cbind(sum(grepl(10, stdy2_dta$post_task2)),"Feel \nFrustration"),
      cbind(sum(grepl(11, stdy2_dta$post_task2)),"Palestinian \nPosts \nare not \nBased \n on Facts"),
      cbind(sum(grepl(12, stdy2_dta$post_task2)),"No \nPatience \nfor \nPalestinians \nPosts")
    )) %>%
  dplyr::rename(.,
                Count = V1,
                Rationale = V2) %>% 
  mutate(., Count = as.numeric(Count))

avoiding <- 
  stdy2_dta %>% 
  filter(.,
         choice_pal!=1)

why_avoid$Share <- as.numeric(why_avoid$Count)/nrow(avoiding)


## Generate histogram----
ggplot(why_avoid, aes(x=reorder(Rationale, -Share), y=Share)) + 
  geom_bar(stat="identity", fill = "dodgerblue4") +
  ylim(0,.6)+
  labs(x = "",
       y = "Share of Respondents")+
  theme_bw()+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 16),
        # axis.text.y = element_text(size=14),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())


#Appendix------
##Figure S5-----
stdy2_dta <-
  stdy2_dta %>% 
  mutate(.,
         arab_therm_bin = ifelse(as.numeric(arab_therm) > 
                                   quantile(as.numeric(stdy2_dta$arab_therm), 0.75),1,0),
         right_wing = ifelse(left_right_id>4,1,0),
         religious = ifelse(rlg>1,1,0))


# Plot percentage of non-avoiders for respondents with different levels of prejudice to Arabs
stdy2_dta %>%
  filter(!is.na(choice_pal_fac)) %>%
  group_by(arab_therm_bin, choice_pal_fac) %>%
  summarise(count = n()) %>%
  mutate(percentage = count / sum(count))  %>% 
ggplot(., aes(x = choice_pal_fac, y = percentage)) +
  geom_bar(stat = "identity", fill = "dodgerblue4") +
  geom_text(aes(label = scales::percent(percentage)), vjust = -0.5, size = 5) +  # Add percentage labels
  labs(x = "Select to Read about Palestinians in Choice Task",
       y = "Percentage of Respondents") +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0, 1)) +  # Show y-axis in percentage format
  facet_grid(~ arab_therm_bin, labeller = as_labeller(c("0" = "Palestinian Thermometer < Top Quartile", 
                                                        "1" = "Palestinian Thermometer Top Quartile"))) +  # Apply custom facet labels
  theme_tufte() +
  theme(text = element_text(size = 16, family = "Times"),
        legend.key = element_blank(),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(size = 12),
        plot.caption = element_text(size = 10, family = "Times", hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        strip.background = element_rect(fill = "lightgray", color = NA),  # Light gray background for facet header
        strip.text = element_text(color = "black", size = 14))  # Customize text color and size


# Plot percentage of non-avoiders by ideology
stdy2_dta %>%
  filter(!is.na(choice_pal_fac)) %>%
  group_by(right_wing, choice_pal_fac) %>%
  summarise(count = n()) %>%
  mutate(percentage = count / sum(count))  %>% 
ggplot(., aes(x = choice_pal_fac, y = percentage)) +
  geom_bar(stat = "identity", fill = "dodgerblue4") +
  geom_text(aes(label = scales::percent(percentage)), vjust = -0.5, size = 5) +  # Add percentage labels
  labs(x = "Select to Read about Palestinians in Choice Task",
       y = "Percentage of Respondents") +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0, 1)) +  # Show y-axis in percentage format
  facet_grid(~ right_wing, labeller = as_labeller(c("0" = "Center & Left Wing", 
                                                    "1" = "Right Wing"))) +  # Apply custom facet labels
  theme_tufte() +
  theme(text = element_text(size = 16, family = "Times"),
        legend.key = element_blank(),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(size = 12),
        plot.caption = element_text(size = 10, family = "Times", hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        strip.background = element_rect(fill = "lightgray", color = NA),  # Light gray background for facet header
        strip.text = element_text(color = "black", size = 14))  # Customize text color and size



# Plot percentage of non-avoiders by religioisity
stdy2_dta %>%
  filter(!is.na(choice_pal_fac)) %>%
  group_by(religious, choice_pal_fac) %>%
  summarise(count = n()) %>%
  mutate(percentage = count / sum(count)) %>% 
ggplot(., aes(x = choice_pal_fac, y = percentage)) +
  geom_bar(stat = "identity", fill = "dodgerblue4") +
  geom_text(aes(label = scales::percent(percentage)), vjust = -0.5, size = 5) +  # Add percentage labels
  labs(x = "Select to Read about Palestinians in Choice Task",
       y = "Percentage of Respondents") +
  scale_y_continuous(labels = scales::percent_format(), limits = c(0, 1)) +  # Show y-axis in percentage format
  facet_grid(~ religious, labeller = as_labeller(c("0" = "Secular", 
                                                   "1" = "Traditional or Religious"))) +  # Apply custom facet labels
  theme_tufte() +
  theme(text = element_text(size = 16, family = "Times"),
        legend.key = element_blank(),
        panel.grid.major = element_blank(), 
        axis.text.x = element_text(size = 12),
        plot.caption = element_text(size = 10, family = "Times", hjust = -.02),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        strip.background = element_rect(fill = "lightgray", color = NA),  # Light gray background for facet header
        strip.text = element_text(color = "black", size = 14))  # Customize text color and size



##Figure S6------
### Create indicators for reasons-----
stdy2_dta <-
  stdy2_dta %>% 
  mutate(.,
         cog_diss = case_when(
           #feelings of discomfort, frustration, despair, fear, or shame
           grepl(1, stdy2_dta$post_task2) ~ 1,
           grepl(10, stdy2_dta$post_task2) ~ 1,
           grepl(4, stdy2_dta$post_task2) ~ 1,
           grepl(9, stdy2_dta$post_task2) ~ 1,
           grepl(3, stdy2_dta$post_task2) ~ 1
         ),
         cog_diss = ifelse(is.na(cog_diss),0, cog_diss),
         naiv_real = case_when(
           #belief in truth or objectivity of one’s own views or belief that outgroup narratives are less truthful or biased; ;
           #anger when exposed to counter-arguments
           grepl(11, stdy2_dta$post_task2) ~ 1,
           grepl(2, stdy2_dta$post_task2) ~ 1 
         ),
         naiv_real = ifelse(is.na(naiv_real),0,naiv_real),
         cog_leth = case_when(
           #no patience; seen enough; other priorities
           grepl(6, stdy2_dta$post_task2) ~ 1,
           grepl(7, stdy2_dta$post_task2) ~ 1,
           grepl(12, stdy2_dta$post_task2) ~ 1 
         ),
         cog_leth = ifelse(is.na(cog_leth),0, cog_leth),
         hopless = case_when(
           #reading posts will not change anything; seen enough
           grepl(5, stdy2_dta$post_task2) ~ 1,
           grepl(6, stdy2_dta$post_task2) ~ 1
         ),
         hopless = ifelse(is.na(hopless),0,hopless))




###Estimate Correlates------
###Cognitive Diss
ideology_cog <-
  lm_robust(cog_diss~ stdize(left_right_id),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Ideology",
         Reason = "Cognitive Dissonance")

religiosity_cog <-
  lm_robust(cog_diss~ stdize(as.numeric(rlg)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Religiosity",
         Reason = "Cognitive Dissonance")


age_group_cog <-
  lm_robust(cog_diss~ stdize(as.numeric(age)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Age Group",
         Reason = "Cognitive Dissonance")


therm_cog <-
  lm_robust(cog_diss~ stdize(as.numeric(arab_therm)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Palastenian \nThermometer",
         Reason = "Cognitive Dissonance")


sm_usage_cog <-
  lm_robust(cog_diss~ stdize(as.numeric(sm_usage)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Frequency on \nSocial Media",
         Reason = "Cognitive Dissonance")


sm_expose_cog <-
  lm_robust(cog_diss~ stdize(as.numeric(sm_exposure)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Value \nDiversity \nOnline",
         Reason = "Cognitive Dissonance")


speak_arabic_cog <-
  lm_robust(cog_diss~ stdize(as.numeric(speak_arabic)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Arabic \nSpeaking",
         Reason = "Cognitive Dissonance")



edu_cog <-
  lm_robust(cog_diss~ stdize(as.numeric(education)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Education",
         Reason = "Cognitive Dissonance")

male_cog <-
  lm_robust(cog_diss~ stdize(as.numeric(male)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Male",
         Reason = "Cognitive Dissonance")


###Naive Realism
ideology_nv <-
  lm_robust(naiv_real~ stdize(left_right_id),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Ideology",
         Reason = "Naive Realism")

religiosity_nv <-
  lm_robust(naiv_real~ stdize(as.numeric(rlg)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Religiosity",
         Reason = "Naive Realism")


age_group_nv <-
  lm_robust(naiv_real~ stdize(as.numeric(age)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Age Group",
         Reason = "Naive Realism")


therm_nv <-
  lm_robust(naiv_real~ stdize(as.numeric(arab_therm)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Palastenian \nThermometer",
         Reason = "Naive Realism")


sm_usage_nv <-
  lm_robust(naiv_real~ stdize(as.numeric(sm_usage)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Frequency on \nSocial Media",
         Reason = "Naive Realism")


sm_expose_nv <-
  lm_robust(naiv_real~ stdize(as.numeric(sm_exposure)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Value \nDiversity \nOnline",
         Reason = "Naive Realism")


speak_arabic_nv <-
  lm_robust(naiv_real~ stdize(as.numeric(speak_arabic)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Arabic \nSpeaking",
         Reason = "Naive Realism")



edu_nv <-
  lm_robust(naiv_real~ stdize(as.numeric(education)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Education",
         Reason = "Naive Realism")

male_nv <-
  lm_robust(naiv_real~ stdize(as.numeric(male)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Male",
         Reason = "Naive Realism")




###Hopeless

ideology_hop <-
  lm_robust(hopless~ stdize(left_right_id),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Ideology",
         Reason = "Hopelessness")

religiosity_hop <-
  lm_robust(hopless~ stdize(as.numeric(rlg)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Religiosity",
         Reason = "Hopelessness")


age_group_hop <-
  lm_robust(hopless~ stdize(as.numeric(age)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Age Group",
         Reason = "Hopelessness")


therm_hop <-
  lm_robust(hopless~ stdize(as.numeric(arab_therm)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Palastenian \nThermometer",
         Reason = "Hopelessness")


sm_usage_hop <-
  lm_robust(hopless~ stdize(as.numeric(sm_usage)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Frequency on \nSocial Media",
         Reason = "Hopelessness")


sm_expose_hop <-
  lm_robust(hopless~ stdize(as.numeric(sm_exposure)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Value \nDiversity \nOnline",
         Reason = "Hopelessness")


speak_arabic_hop <-
  lm_robust(hopless~ stdize(as.numeric(speak_arabic)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Arabic \nSpeaking",
         Reason = "Hopelessness")



edu_hop <-
  lm_robust(hopless~ stdize(as.numeric(education)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Education",
         Reason = "Hopelessness")

male_hop <-
  lm_robust(hopless~ stdize(as.numeric(male)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Male",
         Reason = "Hopelessness")


###Cognitive Leth

ideology_lth <-
  lm_robust(cog_leth~ stdize(left_right_id),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Ideology",
         Reason = "Cognitive Lethargy")

religiosity_lth <-
  lm_robust(cog_leth~ stdize(as.numeric(rlg)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Religiosity",
         Reason = "Cognitive Lethargy")


age_group_lth <-
  lm_robust(cog_leth~ stdize(as.numeric(age)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Age Group",
         Reason = "Cognitive Lethargy")


therm_lth <-
  lm_robust(cog_leth~ stdize(as.numeric(arab_therm)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Palastenian \nThermometer",
         Reason = "Cognitive Lethargy")


sm_usage_lth <-
  lm_robust(cog_leth~ stdize(as.numeric(sm_usage)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Frequency on \nSocial Media",
         Reason = "Cognitive Lethargy")


sm_expose_lth <-
  lm_robust(cog_leth~ stdize(as.numeric(sm_exposure)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Value \nDiversity \nOnline",
         Reason = "Cognitive Lethargy")


speak_arabic_lth <-
  lm_robust(cog_leth~ stdize(as.numeric(speak_arabic)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Arabic \nSpeaking",
         Reason = "Cognitive Lethargy")



edu_lth <-
  lm_robust(cog_leth~ stdize(as.numeric(education)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Education",
         Reason = "Cognitive Lethargy")

male_lth <-
  lm_robust(cog_leth~ stdize(as.numeric(male)),
            data = stdy2_dta) %>% tidy() %>% 
  filter(., 
         term != "(Intercept)") %>% 
  mutate(.,
         term = "Male",
         Reason = "Cognitive Lethargy")

covs_mot <-
  rbind(ideology_cog, religiosity_cog, age_group_cog, therm_cog,sm_usage_cog, sm_expose_cog, speak_arabic_cog, edu_cog, male_cog,
        ideology_nv, religiosity_nv, age_group_nv, therm_nv, sm_usage_nv, sm_expose_nv, speak_arabic_nv, edu_nv, male_nv,
        ideology_hop, religiosity_hop, age_group_hop, therm_hop,sm_usage_hop, sm_expose_hop, speak_arabic_hop, edu_hop, male_hop,
        ideology_lth, religiosity_lth, age_group_lth, therm_lth, sm_usage_lth, sm_expose_lth, speak_arabic_lth, edu_lth, male_lth)


ggplot(covs_mot, aes(x=estimate, y=reorder(term, estimate))) + 
  geom_point(color = "dodgerblue4")+
  geom_errorbar(aes(xmin=conf.low, xmax=conf.high), 
                width=.1, color = "dodgerblue4") +
  geom_vline(xintercept = 0,
             color = "black",
             linetype = "dashed",
             size = 0.4) +
  facet_wrap(~Reason)+
  xlim(-.15,.15)+
  labs(x = "Bivariate Correlation",
       y = "")+
  scale_color_manual(values = c("dodgerblue4"))+
  theme_bw()+
  theme(panel.spacing = unit(1.5, "lines"),
        text = element_text(size = 14),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank(),
        legend.key=element_blank())


#Table S9------
# Prepare data for descriptive stats table
covs <- 
  stdy2_dta %>%
  mutate(.,
         Age = case_when(
           age == 1 ~ "18-24",
           age == 2 ~ "25-34",
           age == 3 ~ "35-44",
           age == 4 ~ "45-54",
           age == 5 ~ "55-64",
           age == 6 ~ "65-70"
         ),
         Gender = case_when(
           male == 1 ~ "Male",
           female == 1 ~ "Female"
         ),
         Religiosity = case_when(
           rlg == 1 ~ "Secular",
           rlg == 2 ~ "Traditional",
           rlg == 3 ~ "Religious",
           rlg == 4 ~ "Ultra-Orthodox"
         )) %>% 
  dplyr::select(Age, Gender, Religiosity)

# Calculate proportions excluding "NAs"
proportions_df <- covs %>%
  pivot_longer(cols = everything(), names_to = "Variable", values_to = "Value") %>%
  filter(!is.na(Value)) %>%  # Exclude "No-Response"
  group_by(Variable, Value) %>%
  summarise(
    Count = n(),  # Count of each value
    .groups = "drop_last"  # Retain grouping by Variable for next step
  ) %>%
  mutate(
    Total = sum(Count),  # Total counts for each variable
    Proportion = Count / Total  # Proportion for each value
  ) %>%
  ungroup() %>%  # Remove all grouping
  arrange(Variable, desc(Proportion)) %>% 
  dplyr:: select(Variable,Value,Proportion) %>% 
  mutate(.,
         Proportion = round(Proportion,3))

# read in ipanel data
ipanel <- read.csv("ipanel_quota.csv")
`Ipanel prop.` <-ipanel$IPanel.Proportion
descriptive_table <- cbind(proportions_df,`Ipanel prop.`)
  
# Save the LaTeX table
descriptive_table %>%
  kable(format = "latex", 
        booktabs = TRUE, 
        align = "c", 
        label = "tab:desc_2",
        caption = "Descriptive Statistics Study 2") 


